// -*- coding: utf-8 -*-

/**
 * @file hybrid_ode_solver.h
 * @brief A wrapper for HODES.
 */

#ifndef _HYBRID_ODE_SOLVER_H
#define _HYBRID_ODE_SOLVER_H

#define PIPE_CACHING

#include <iostream>
#include <memory>
#include <set>
#include <map>

#ifdef USE_PIPE_HASHMAP
#ifndef __GNUC__
#include <hash_map>
#else
#include <ext/hash_map>
#endif
#endif

#include "vnode.h"

#include "hybrid_config.h"
#include "hybrid_var.h"
#include "hybrid_model.h"
#include "hybrid_util.h"

NAMESPACE_ELISA

/// the minimum step size for solving ODEs.
extern const double MinH;
/// the ratio for MinH versus the accracy of the time domain approximation.
extern const int    Ratio;

/**
 * @ingroup Hybrid
 * @class HybridOdeSolver
 */
class HybridOdeSolver
{
public:
	/*
	HybridOdeSolver(vnodelp::FADBAD_AD *ad, VNInterval t0, VNIntVec x0)
		: m_odeSolver(ad),
		  m_initTime(t0),
		  m_initValue(x0)
	{
		setup();
	}
	*/
	
	HybridOdeSolver(ContConstraint & cc, bool dump = false)
		: m_odeSolver(dynamic_cast<ContConstraintRep *>(cc.rep())->getAD().get()),
		  m_constr(cc)
	{
		HybridOdeSolver::dump = dump;
		if (dump) {
			odeOut << std::setprecision(20);
			//startDumpTightEncl();
			//atexit(endDumpTightEncl);
		}

		// 制約中の時間変数にアクセスする
		ContConstraintRep *rep = cc.rep();
		Time& time = rep->getTime();

		// 初期状態の集合
		//m_initValue = time.getInitValue(); // これでよいはずでは。。
		for (VNIntVec::iterator it(time.getInitValue().begin());
			 it != time.getInitValue().end(); ++it) {

			m_initValue.push_back(*it);
		}		
		// 初期時間の集合
		m_initTime  = time.getInitTime();

		setup();
	}

	/// ある時点における状態の包囲を返す。
	virtual const VNIntVec getEnclosure(double t) = 0;
	/// ある時間区間における状態の包囲を返す。
	virtual const VNIntVec getEnclosure(const VNInterval& t) = 0;
	/// ある時間区間の端でとった最小時間区間における状態の包囲を返す。
	virtual void getBBox(bool left, const VNInterval& t, VNInterval& bT, VNIntVec& bAX);
	/// ある時間区間の中心でとった最小時間区間における状態の包囲を返す。
	virtual void getMEncl(const VNInterval& t, double& midT, VNIntVec& midAX);

	const VNInterval& getInitTime() {
		return m_initTime;
	}
	const VNIntVec& getInitValue() {
		return m_initValue;
	}

	vnodelp::VNODE& getVnode() {
		return m_odeSolver;
	}

	virtual void savePipe(VNPipe pipe) = 0;

	void setupLO() {
		setup();
		m_odeSolver.setOrder(ODE_ORDER_LO);
	}

	static bool dump;

protected:
	void setup() {
		m_odeSolver.setTols(ODE_ATOL, ODE_RTOL);
		m_odeSolver.setOrder(ODE_ORDER);
		m_odeSolver.setHmin(MinH);
		//m_odeSolver.setOneStep(vnodelp::on);
	}

	vnodelp::VNODE m_odeSolver;

	VNInterval	m_initTime;
	VNIntVec	m_initValue;

	ContConstraint& m_constr;
};

class HybridOdeSolverImpl : public HybridOdeSolver
{
	struct PipeCmp
	{
		/// 比較関数は, 等しい値に対して常にfalseを返すようにする (厳格弱順序)
		bool operator()(VNIntPtr l, VNIntPtr r) const {
			if (v_bias::sup(*l) < v_bias::sup(*r))
				return true;
			else
				if (v_bias::sup(*l) == v_bias::sup(*r))
					return (v_bias::width(*l) > v_bias::width(*r));
				else
					return false;
		}

		/*
		bool operator()(const VNIntPtr l, const VNIntPtr r) const {
		if (v_bias::sup(*l) < v_bias::sup(*r))
		return true;
		else
		if (v_bias::sup(*l) == v_bias::sup(*r))
		return (v_bias::width(*l) < v_bias::width(*r));
		}
		*/
	};

	typedef std::map<VNIntPtr, VNEncl, PipeCmp> PipeSet;

public:
	/*
	HybridOdeSolverImpl(vnodelp::FADBAD_AD *ad, VNInterval t0, VNIntVec x0)
		: HybridOdeSolver(ad, t0, x0)
	{
		initPipeSet();
	}
	*/

	HybridOdeSolverImpl(ContConstraint & cc, bool dump = false)
		: HybridOdeSolver(cc, dump)
	{
		initPipeSet();
	}

	/// ある時点における状態の包囲を返す。
	virtual const VNIntVec getEnclosure(double t);
	/// ある時間区間における状態の包囲を返す。
	virtual const VNIntVec getEnclosure(const VNInterval& t);

	virtual void savePipe(VNPipe pipe) {
		m_pipeSet.insert(pipe);
	}

private:
	void initPipeSet() {
		if (v_bias::width(m_initTime) == 0) {
			VNIVPtr		x0p(new VNIntVec(m_initValue));
			VNIntPtr	t0p(new VNInterval(m_initTime));
			m_pipeSet.insert(VNPipe(t0p, VNEncl(x0p, x0p)));
		}
	}

	/// 時間区間 (上端値) をキーにしたパイプの集合。
	PipeSet		m_pipeSet;
};


#ifdef USE_PIPE_HASHMAP
class HybridOdeSolverImpl1 : public HybridOdeSolver
{
	/*
	struct PipeHF
	{
	bool operator()(VNIntPtr k) const {
	return std::hash<double>(v_bias::sup(k));
	}
	};

	struct PipeEK
	{
	bool operator()(VNIntPtr k) const {
	return std::hash<double>(v_bias::sup(k));
	}		
	};
	*/

#ifndef __GNUC__
	typedef stdext::hash_multimap<double, VNPipe>		PipeSet;
#else
	struct hashString { // the name is misleading - a changement of other code 
		inline std::size_t operator()(const double& key) const {
			return ((size_t)key);
		}
	};
	typedef __gnu_cxx::hash_multimap<double, VNPipe, hashString> PipeSet;
#endif
	typedef std::pair<double, VNPipe>			PipeExt;

public:
	/*
	HybridOdeSolverImpl1(vnodelp::FADBAD_AD *ad, VNInterval t0, VNIntVec x0)
		: HybridOdeSolver(ad, t0, x0)
		  //, m_pipeSet(1024)
	{
		initPipeSet();
	}
	*/

	HybridOdeSolverImpl1(ContConstraint & cc, bool dump = false)
		: HybridOdeSolver(cc, dump)
		  //, m_pipeSet(1024)
	{
		initPipeSet();
	}

	/// ある時点における状態の包囲を返す。
	virtual const VNIntVec getEnclosure(double t);
	/// ある時間区間における状態の包囲を返す。
	virtual const VNIntVec getEnclosure(const VNInterval& t);

	virtual void savePipe(VNPipe pipe) {
		m_pipeSet.insert(PipeExt(v_bias::sup(*pipe.first), pipe));
		//m_pipeSet[v_bias::sup(*pipe.first)] = pipe;
	}

private:
	void initPipeSet() {
		if (v_bias::width(m_initTime) == 0) {
			VNIVPtr		x0p(new VNIntVec(m_initValue));
			VNIntPtr	t0p(new VNInterval(m_initTime));

			m_pipeSet.insert(PipeExt(v_bias::sup(m_initTime), VNPipe(t0p, VNEncl(x0p, x0p))));
			//VNPipe pipe(t0p, VNEncl(x0p, x0p));
			//m_pipeSet[1.0] = pipe;
		}
	}

	/// 時間区間 (上端値) をキーにしたパイプの集合。
	PipeSet		m_pipeSet;
};
#endif // USE_PIPE_HASHMAP

END_NAMESPACE_ELISA

#endif // HYBRID_ODE_SOLVER
